function queue_dynamics_delayed(provisioning_time,hold_duration)
randn('seed',12345)
lambda = randn(1,400);
%lambda(100:200)=4+lambda(100:200);
%lambda=lambda+10;
% randn('seed',16576576576);  k=cumsum(randn(200,1)*5 ); k= k-min(k); plot(k)
lambda = lambda + sin((1:400)/20)*10+30;
lambda = [0 lambda];
lambda(181:end)=0;
a = 1;
b = [1/4 1/4 1/4 1/4];
lambda_smooth = filter(b,a,lambda);


r=1; target_q=10;
x0=20; 


% prices_amz().get_hourly_price()./prices_amz().get_compute_units()
inst0 = [0;0;0];


%provisioning_time=5; %5
%hold_duration=1;   %60
% A=[0 0 0
%       1 0 0
%       0 1 0];
stations=provisioning_time+hold_duration;
A=[zeros(1,stations-1)
      eye(stations-1)];
A=[A zeros(stations,1)];   

% Ain = [1
%             0
%             0 ]; 
Ain=zeros(stations,1); 
Ain(1,1)=1; 

 T=10;
% %m=mean(lambda(:,2:T+1));
 Vartheta0=zeros(stations,1); 
 Vartheta0(provisioning_time+1,1) = 5; 
% cvx_begin
%     variables X(1,T+1) Vartheta(stations,T+1) U(1,T)  cost_r
%     Vartheta(:,2:T+1) == A*Vartheta(:,1:T)+Ain*U(:,1:T);
%     X(:,2:T+1) == X(:,1:T)+lambda(:,2:T+1)-[zeros(1,provisioning_time) ones(1,hold_duration)]*Vartheta(:,1:T)  ; 
%     X(:,1) == x0;      
%     Vartheta(:,1)== Vartheta0;  % ones(stations,1); %[10; 4; 4]; 
%     Vartheta>0; 
%      X(:,2:T+1)>0;
%     cost_r==r*sum(U(:,1:T));
%     minimize (sum_square(target_q*ones(1,T)-X(:,2:T+1))...
%                         +cost_r)
% cvx_end
% Xopt=X;
% Uopt=U; 
% optcost = sum(sum_square(Xopt(:,2:T+1)-target_q))...
%                         +r*sum(Uopt)
% 
% subplot(3,1,1); stairs(Uopt); title('Uopt'); 
% subplot(3,1,2); plot(Xopt); title('Q opt');
% subplot(3,1,3); plot(lambda_smooth); title('lambda'); 
% figure; plot(([zeros(1,provisioning_time) ones(1,hold_duration)]*Vartheta)'); 
% subplot(4,1,1); plot(Xopt); title('Xopt') 
% subplot(4,1,2); plot(Uopt'); title('Uopt') 
% subplot(4,1,3); plot(Xall); title('Xall') 
% subplot(4,1,4); plot(Uall'); title('Uall') 
%
% subplot(3,1,3); plot(lambda(2:end))

kk=0; 

%clear
% MPC
for T=1:60 %
%T = 60; % horizon
nsteps = 180; % number of steps 
Vartheta0=zeros(stations,1); 
Vartheta0(provisioning_time+1,1) = 5;  
n=1; m=1; 
x = x0;  vartheta=Vartheta0;

Xall = zeros(n,nsteps); Uall = zeros(m,nsteps); VarthetaAll=zeros(1,nsteps); 
alpha_est= 0.75;
v_cost=10; 
lambda_hat=0.1576;
D=[zeros(1,provisioning_time) ones(1,hold_duration)];
dur = [];
for i = 1:nsteps
    fprintf('.');  
    tic
    cvx_begin quiet 
        variables X(1,T+1) U(1,T)  Vartheta(stations,T+1)  cost_r  lmbda(1,T)
        % here I substituted lambda(:,2:T+1) with [mean(lambda(100:200)) mean(lambda(100:200))]
        X(:,1) == x;           
       lmbda == lambda(:,1+i:T+i); %lambda_hat*ones(1,T)
        X(:,2:T+1) == X(:,1:T)+ lmbda -D*Vartheta(:,1:T);        
         Vartheta(:,1)== vartheta;  % ones(stations,1); %[10; 4; 4]; 
         Vartheta(:,2:T+1) == A*Vartheta(:,1:T)+Ain*U(:,1:T);
         Vartheta>0; 
         %X(:,2:T+1)>0;
         cost_r==r*sum(U(:,1:T));
        minimize (sum_square(target_q*ones(1,T)-X(:,2:T+1))...
                         +cost_r...
                         +v_cost *  sum_square(X(:,T+1)-target_q))
    cvx_end
    dur(i) = toc;
    Xall(:,i) = x; u = U(:,1); Uall(:,i) = u;
    x = x+ lambda(:,i) -D*vartheta;   
    vartheta = A*vartheta+Ain*u;
    VarthetaAll(:,i) = D*vartheta; 
    lambda_hat= alpha_est*lambda_hat+(1-alpha_est)*lambda(:,i+1);
end

mpccost = sum(sum_square(Xall-target_q))...
                        +r*sum(Uall);
% disp(sprintf('T=%d',T)); disp(dur); 
disp(sprintf('mpccost=%d',mpccost));               
disp(sprintf('cvx_slvitr=%d',cvx_slvitr));
end
%  plot(Uall)
% figure
% plot(Xall);

kk=0

exit





